## This file replicates the figures and tables in Chapter 12 of:
## Sharan Grewal, Soldiers of Democracy? Military Legacies and
## the Arab Spring, Oxford University Press, 2023.
## For any questions, contact: ssgrewal@wm.edu



######################
## Load Survey Data ##
######################

survey1 <- read.csv("survey1.csv")
survey2 <- read.csv("survey2.csv")
survey3 <- read.csv("survey3.csv")

## Load packages
library(ggplot2)
library(scales)
library(stargazer)
library(effects)
library(sjPlot)
library(list)
library(xtable)



#############
## Figures ##
#############

## Survey 1-Rank

survey1$rankTUN <- factor(survey1$rankTUN, levels=(c("Private", "Private First Class","Corporal", "Master Corporal", "Sergeant", "Staff Sergeant", "Sergeant First Class","Master Sergeant", "Sergeant Major", "Second Lieutenant", "First Lieutenant", "Captain", "Major", "Lieutenant Colonel", "Colonel", "Colonel Major","Brigadier General", "Divisional General", "Army Corps General")))

par(mar=c(10,4,4,2))
barplot <- barplot((table(survey1$rankTUN)/nrow(survey1[!is.na(survey1$rankTUN),]))*100, ylab="Percent", xlab="", main="Survey 1: Military Ranks, Tunisia (N=72)", ylim=c(0,70),las=2)
text(x=t(barplot), y=(table(survey1$rankTUN)/nrow(survey1[!is.na(survey1$rankTUN),]))*100, labels=(round((table(survey1$rankTUN)/nrow(survey1[!is.na(survey1$rankTUN),]))*100,0)), pos=3)


## Survey 2-Rank

survey2$rankTUN <- factor(survey2$rankTUN, levels=(c("Private", "Private First Class","Corporal", "Master Corporal", "Sergeant", "Staff Sergeant", "Sergeant First Class","Master Sergeant", "Sergeant Major", "Second Lieutenant", "First Lieutenant", "Captain", "Major", "Lieutenant Colonel", "Colonel", "Colonel Major","Brigadier General", "Divisional General", "Army Corps General")))
survey2$rank3 <- factor(survey2$rank3, levels=(c("Anonymous","Soldier","NCO","Jun. Officer","Sen. Officer")))

par(mar=c(10,4,4,2))
barplot <- barplot((table(survey2$rankTUN)/nrow(survey2[!is.na(survey2$rankTUN),]))*100, ylab="Percent", xlab="", main="Survey 2: Military Ranks, Tunisia (N=194)", ylim=c(0,50),las=2)
text(x=t(barplot), y=(table(survey2$rankTUN)/nrow(survey2[!is.na(survey2$rankTUN),]))*100, labels=(round((table(survey2$rankTUN)/nrow(survey2[!is.na(survey2$rankTUN),]))*100,0)), pos=3)


## Rank-Egypt

survey3$rankEG <- factor(survey3$rankEG, levels=(c("Private", "Corporal", "Sergeant", "Staff Sergeant", "Warrant Officer", "Chief Warrant Officer", "Second Lieutenant", "Lieutenant", "Captain", "Major", "Lieutenant Colonel", "Colonel", "Brigadier", "Major General", "Lieutenant General", "General", "Field Marshall")))
survey3$rank3 <- factor(survey3$rank3, levels=(c("Anonymous","Soldier","NCO","Jun. Officer","Sen. Officer")))

par(mar=c(10,4,4,2))
barplot <- barplot((table(survey3$rankEG)/nrow(survey3[!is.na(survey3$rankEG),]))*100, ylab="Percent", xlab="", main="Survey 3: Military Ranks, Egypt (N=1668)", ylim=c(0,50),las=2)
text(x=t(barplot), y=(table(survey3$rankEG)/nrow(survey3[!is.na(survey3$rankEG),]))*100, labels=(round((table(survey3$rankEG)/nrow(survey3[!is.na(survey3$rankEG),]))*100,0)), pos=3)


dev.off()




####################
## Ind. Variables ##
####################


## Survey 1
data.summary <- cbind(rbind(mean(survey1$neglect<3, na.rm=T),
mean(survey1$neglect>3, na.rm=T)),
rbind(mean(survey1$favorMOI<3, na.rm=T),
mean(survey1$favorMOI>3, na.rm=T)),
rbind(mean(survey1$sahel<3, na.rm=T),
mean(survey1$sahel>3, na.rm=T)))*100

barplot <- barplot(data.summary, ylim=c(0,100), beside=T, ylab="Percent", main="Survey of Tunisian Officers' Association (N=72)", xaxt = "n", legend.text=c("Disagree","Agree"), las=2)
text(x=t(barplot), y=t(data.summary), labels=(round(t(data.summary),0)), pos=3)
axis(1, at=c(2,5,8), labels=c("Ben Ali\nneglected\nthe military", "Ben Ali\nfavored\nthe MOI", "Ben Ali\nfavored\nthe Sahel"), tick=F, line=2)


## Survey 2
data.summary <- cbind(rbind(mean(survey2$BA_mil_1<3, na.rm=T),
mean(survey2$BA_mil_1>3, na.rm=T)),
rbind(mean(survey2$BA_mil_2<3, na.rm=T),
mean(survey2$BA_mil_2>3, na.rm=T)),
rbind(mean(survey2$BA_mil_3<3, na.rm=T),
mean(survey2$BA_mil_3>3, na.rm=T)))*100

barplot <- barplot(data.summary, ylim=c(0,100), beside=T, ylab="Percent", main="Online Survey of Tunisian Military (N=271)", xaxt = "n", legend.text=c("Disagree","Agree"), las=2, args.legend=list(y=103))
text(x=t(barplot), y=t(data.summary), labels=(round(t(data.summary),0)), pos=3)
axis(1, at=c(2,5,8), labels=c("Ben Ali\nneglected\nthe military", "Ben Ali\nfavored\nthe MOI", "Ben Ali\nfavored\nthe Sahel"), tick=F, line=2)


## Survey 3
data.summary <- cbind(rbind(mean(survey3$mub_mil_1<3, na.rm=T),
mean(survey3$mub_mil_1>3, na.rm=T)),
rbind(mean(survey3$mub_mil_2<3, na.rm=T),
mean(survey3$mub_mil_2>3, na.rm=T)),
rbind(mean(survey3$mub_mil_3<3, na.rm=T),
mean(survey3$mub_mil_3>3, na.rm=T)))*100

barplot <- barplot(data.summary, ylim=c(0,100), beside=T, ylab="Percent", main="Online Survey of Egyptian Military (N=2171)", xaxt = "n", legend.text=c("Disagree","Agree"), las=2, args.legend=list(y=103))
text(x=t(barplot), y=t(data.summary), labels=(round(t(data.summary),0)), pos=3)
axis(1, at=c(2,5,8), labels=c("Mubarak\nneglected\nthe military", "Mubarak\nfavored\nthe MOI", "Gamal\nwas a\nthreat"), tick=F, line=2)



################
## Transition ##
################

## Survey 1: Interests
table(survey1$int)


## Survey 2: Troika
data.summary <- cbind(rbind(mean(survey2$troika_mil_2<3, na.rm=T),
mean(survey2$troika_mil_2>3, na.rm=T)),
rbind(mean(survey2$troika_mil_1<3, na.rm=T),
mean(survey2$troika_mil_1>3, na.rm=T)),
rbind(mean(survey2$troika_mil_3<3, na.rm=T),
mean(survey2$troika_mil_3>3, na.rm=T)))*100

barplot <- barplot(data.summary, ylim=c(0,100), beside=T, ylab="Percent", main="Online Survey of Tunisian Military (N=271)", xaxt = "n", legend.text=c("Disagree","Agree"), las=2, args.legend=list(y=103))
text(x=t(barplot), y=t(data.summary), labels=(round(t(data.summary),0)), pos=3)
axis(1, at=c(2,5,8), labels=c("Troika\nenhanced\nmaterially", "Troika\nenhanced\npolitically", "Troika\nended\ndiscrimination"), tick=F, line=2)


## Survey 3: Morsi
data.summary <- cbind(rbind(mean(survey3$morsi_mil_1<3, na.rm=T),
mean(survey3$morsi_mil_1>3, na.rm=T)),
rbind(mean(survey3$morsi_mil_2<3, na.rm=T),
mean(survey3$morsi_mil_2>3, na.rm=T)),
rbind(mean(survey3$morsi_mil_3<3, na.rm=T),
mean(survey3$morsi_mil_3>3, na.rm=T)),
rbind(mean(survey3$morsi_mil_4<3, na.rm=T),
mean(survey3$morsi_mil_4>3, na.rm=T)))*100

barplot <- barplot(data.summary, ylim=c(0,100), beside=T, ylab="Percent", main="Online Survey of Egyptian Military (N=2171)", xaxt = "n", legend.text=c("Disagree","Agree"), las=2, args.legend=list(y=103))
text(x=t(barplot), y=t(data.summary), labels=(round(t(data.summary),0)), pos=3)
axis(1, at=c(2,5,8,11), labels=c("Morsi\nbad security\ndecisions", "Morsi\ndenied\ncontracts","Morsi\nreduced\nconst. powers","Morsi\nBrotherhoodized\nmilitary"), tick=F, line=2)




#################
## Composition ##
#################

## Class

t.test(survey2$class[survey2$rank2=="Officer"]==5, survey3$class[survey3$officer==1]==5)

data.summary <- rbind(table(survey2$class[survey2$officer==1])/59,
table(survey3$class[survey3$officer==1])/275)*100

barplot <- barplot(data.summary,
main="What Class are Military Officers? (among Officers)", xaxt="n", cex.names=0.9, ylab="Percent", ylim=c(0, 62), beside=T, legend.text=c("Tunisia","Egypt"), args.legend=list(x=3.5, y=62))
axis(1, at=c(2,5,8,11,14), labels=c("Lower\nclass","Lower-middle\nclass","Middle\nclass","Upper-middle\nclass","Upper\nclass"), tick=F, line=1)
text(barplot, data.summary, labels=round(data.summary, 0), pos=3)
text(c(5,8,11,14), 53, labels=c("***","***","***","***"), pos=3, cex=1.5)



## Why join?

t.test(survey2$join_wealth[survey2$rank2=="Officer"], survey3$join_wealth[survey3$rank2=="Officer"])

data.summary <- cbind(rbind(mean(survey2$join_com[survey2$officer==1], na.rm=T),
      mean(survey3$join_com[survey3$officer==1], na.rm=T)),
  rbind(mean(survey2$join_pri[survey2$officer==1], na.rm=T),
      mean(survey3$join_pri[survey3$officer==1], na.rm=T)),
  rbind(mean(survey2$join_mob[survey2$officer==1], na.rm=T),
      mean(survey3$join_mob[survey3$officer==1], na.rm=T)),
  rbind(mean(survey2$join_wealth[survey2$officer==1], na.rm=T),
      mean(survey3$join_wealth[survey3$officer==1], na.rm=T)))*100

barplot <- barplot(data.summary,
main="Reasons for Joining the Military (among Officers)", xaxt="n", cex.names=0.9, ylab="Percent", ylim=c(0, 60), beside=T, legend.text=c("Tunisia","Egypt"), args.legend=list(x=3, y=60))
axis(1, at=c(2,5,8,11), labels=c("Work in military\ncompany after","Pays more than\nprivate sector","Upward\nmobility","Total\n"), tick=F, line=1)
text(barplot, data.summary, labels=round(data.summary, 0), pos=3)
text(c(2,5,8,11), 35, labels=c("**","***","**","***"), pos=3, cex=1.5)





#####################
## Professionalism ##
#####################


## how define professionalism

t.test(survey2$profdef_2[survey2$officer==1]==1, survey3$profdef_2[survey3$officer==1]==1)

data.summary <- cbind(rbind(mean(survey2$profdef_2[survey2$officer==1]==1, na.rm=T), mean(survey3$profdef_2[survey3$officer==1]==1, na.rm=T)),
rbind(mean(survey2$profdef_1[survey2$officer==1]==1, na.rm=T), mean(survey3$profdef_1[survey3$officer==1]==1, na.rm=T)),
rbind(mean(survey2$profdef_3[survey2$officer==1]==1, na.rm=T), mean(survey3$profdef_3[survey3$officer==1]==1, na.rm=T)),
rbind(mean(survey2$profdef_6[survey2$officer==1]==1, na.rm=T), mean(survey3$profdef_6[survey3$officer==1]==1, na.rm=T)),
rbind(mean(survey2$profdef_4[survey2$officer==1]==1, na.rm=T), mean(survey3$profdef_4[survey3$officer==1]==1, na.rm=T)),
rbind(mean(survey2$profdef_5[survey2$officer==1]==1, na.rm=T), mean(survey3$profdef_5[survey3$officer==1]==1, na.rm=T)))*100
barplot <- barplot(data.summary,
main="What Does it Mean to be Professional? (among Officers)",
names.arg=c("Being\nApolitical", "Technical\nExpertise", "Military int.\nabove own", "Meritocratic\nAppointments", "Protecting\nNat Sec","Following\nOrders"), cex.names=0.82, ylab="Most important factor (%)", ylim=c(0, 60), beside=T,legend.text=c("Tunisia","Egypt"), xpd=F)
text(barplot, data.summary, labels=round(data.summary, 0), pos=3)
text(2, 40, labels=c("**"), pos=3, cex=1.5)





####################
## Dep. Variables ##
####################


## DVs: Association

data.summary <- cbind(rbind(mean(survey1$benali<3, na.rm=T),
mean(survey1$benali==3, na.rm=T), mean(survey1$benali>3, na.rm=T)),
rbind(mean(survey1$democracy<3, na.rm=T), mean(survey1$democracy==3, na.rm=T),
mean(survey1$democracy>3, na.rm=T)))*100

barplot <- barplot(data.summary, ylim=c(0,100), beside=T, ylab="Percent", main="Survey of Tunisian Officers' Association (N=72)", xaxt = "n", legend.text=c("Oppose","Neutral","Support"), las=2, args.legend=list(y=103, x=5))
text(x=t(barplot), y=t(data.summary), labels=(round(t(data.summary),0)), pos=3)
axis(1, at=c(2.5, 6.5), labels=c("Attitudes toward Ben Ali", "Attitudes toward Democracy"), tick=F, line=1)



## DV: Fire

t.test(survey2$fire_BA[survey2$officer==1], survey3$fire[survey3$officer==1])

data.summary <- cbind(rbind(mean(survey2$fire_BA, na.rm=T), mean(survey3$fire, na.rm=T)),
rbind(mean(survey2$fire_BA[survey2$officer!=1], na.rm=T), mean(survey3$fire[survey3$officer!=1], na.rm=T)),
rbind(mean(survey2$fire_BA[survey2$officer==1], na.rm=T),
mean(survey3$fire[survey3$officer==1], na.rm=T)))*100

barplot <- barplot(data.summary, beside=T, ylim=c(0,100), ylab="Percent", main="Would Have Fired on Protesters in 2011", xaxt = "n", las=2)
text(x=t(barplot), y=t(data.summary), labels=(round(t(data.summary),0)), pos=3)
text(x=c(2,5,8), y=c(40,40,40), labels=c("","","*"), pos=3, cex=1.5)
axis(1, at=c(1.5,2.5,4.5,5.5,7.5,8.5), labels=c("Tunisia\nOnline\n(N=271)", "Egypt\nOnline\n(N=2171)","Tunisia\nSoldiers\n(N=209)", "Egypt\nSoldiers\n(N=1888)","Tunisia\nOfficers\n(N=62)","Egypt\nOfficers\n(N=283)"), tick=F, line=1, cex.axis=0.79)



## DV: Coup

data.summary <- cbind(rbind(mean(survey2$coup_ammar>3, na.rm=T), mean(survey3$direct, na.rm=T)),
rbind(mean(survey2$coup_ammar[survey2$officer!=1]>3, na.rm=T), mean(survey3$direct[survey3$officer!=1], na.rm=T)),
rbind(mean(survey2$coup_ammar[survey2$officer==1]>3, na.rm=T),
mean(survey3$direct[survey3$officer==1], na.rm=T)))*100

barplot <- barplot(data.summary, beside=T, ylim=c(0,100), ylab="Percent", main="Support for Military Coup in 2013", xaxt = "n", las=2)
text(x=t(barplot), y=t(data.summary), labels=(round(t(data.summary),0)), pos=3)
text(x=c(2,5,8), y=c(90,90,90), labels=c("***","***","***"), pos=3, cex=1.5)
axis(1, at=c(1.5,2.5,4.5,5.5,7.5,8.5), labels=c("Tunisia\nOnline\n(N=271)", "Egypt\nOnline\n(N=2171)","Tunisia\nSoldiers\n(N=209)", "Egypt\nSoldiers\n(N=1888)","Tunisia\nOfficers\n(N=62)","Egypt\nOfficers\n(N=283)"), tick=F, line=1, cex.axis=0.79)




#################
## Regressions ##
#################


## Survey 1: Ben Ali and Democracy

one <- lm(benali~counterbalanced+favor_sahel+fromSahel+islamist+train_US+army+(rank3), data=survey1)
two <- lm(democracy~prime+improve+favor_sahel+fromSahel+islamist+train_US+army+(rank3), data=survey1)

stargazer(one, two, 
      dep.var.labels=c("Support for Ben Ali","Support for Democracy"), title="Survey of Tunisian Officers' Association",
      covariate.labels=c("Ben Ali counterbalanced","Prime","Troika rebalanced","Ben Ali favored Sahel","from Sahel","Islamist","US Training","Army","Rank"))


## Plots

one <- lm(benali~counterbalance+sahel+fromSahel+islamist+train_US+army+(rank3), data=survey1)
two <- lm(democracy~prime+int+sahel+fromSahel+islamist+train_US+army+(rank3), data=survey1)

plot(effect("counterbalance", one), rug=F, xlab="Ben Ali Counterbalanced the Military (1-5)", ylab="Support for Ben Ali (1-5)", main="Counterbalancing and Resentment toward Ben Ali", colors="black")

plot(effect("int", two), rug=F, xlab="Troika Rebalanced the Military (1-5)", ylab="Support for Democracy (1-5)", main="Rebalancing and Support for Democracy", colors="black")

five <- effect("fromSahel", one, xlevels=list(fromSahel=c(0,1)))
data.summary <- data.frame("fit"=five$fit, "se.fit"=five$se, "name"=c("Interior","Coast"))

ggplot(data.summary, 
       aes(x=name, y=fit)) +  
  geom_bar(position=position_dodge(width=0.5), stat="identity", fill=c("grey40","grey80")) + 
  geom_errorbar(aes(ymin=fit-1.96*se.fit, ymax=fit+1.96*se.fit, width=0.1), size=0.5) +
  ggtitle("Regional Origins and Support for Ben Ali") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(1,4),breaks=seq(1,4,by=0.5),oob=rescale_none) +
  xlab("") +
  ylab("Support for Ben Ali") +
  geom_text(data=data.summary, 
            aes(x=name, y=fit,
            label=formatC(round(fit,2),format='f',digits=2)),
            vjust=-3,size=5.5,position=position_dodge(0.9)) +
  geom_text(aes(x=1.5, y=3.8),size=5,label="p=0.014") +
  theme(text = element_text(size=17))


plot(effect("sahel", two), rug=F, xlab="Ben Ali favored Sahel (1-5)", ylab="Support for Democracy (1-5)", main="Regional Favoritism and Democracy", colors="black")




##########
## Fire ##
##########

one <- lm((fire_BA)~counterbalanced+Region+eliteoff+join_wealth+join_gov+apol+cons+islamist+pray4+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey2)

two <- lm((fire_BA)~favor_sahel*Region+eliteoff+join_wealth+join_gov+apol+cons+islamist+pray4+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey2)

three <- lm((fire)~counterbalanced+eliteoff+join_wealth+join_gov+apol+cons+islamist+pray4+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey3)

four <- lm((fire)~gamal+eliteoff+join_wealth+join_gov+apol+cons+islamist+pray4+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey3)

stargazer(one, two, three, four, single.row=T, digits=2, covariate.labels=c("Counterbalanced","Ben Ali favored coast","from coast","Gamal was threat","Elite Officer","Joined for Wealth","Joined for Power","Apolitical","Conscript","Islamist","Prayer","Registered to Vote","Western Training","Joined to Defend Country","Joined for Thrill of Fighting","Deployed 2011","Active-Duty","Army","Navy","Soldier","NCO","Jun. Officer","Sen. Officer","Education","Rural","Ben Ali favored coast*from coast"), dep.var.labels=c("Ben Ali (Tunisia)","Mubarak (Egypt)"), dep.var.caption="Dependent variable: Would have fired to protect...", title="Willingness to Fire on Protesters in 2011 Revolution")


## Plots 

one <- lm((fire_BA)~counterbalance+Region+eliteoff+join_wealth+join_gov+apol+cons+islamist+pray4+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey2)
plot(effect("counterbalance", one), main="Willingness to fire on 2011 protesters (Tunisia)", xlab="Ben Ali counterbalanced the military (5=str. agree)", ylab="Would have fired (0-1)", colors="black", band.colors="grey0")

two <- lm((fire_BA)~BA_mil_3*Region+eliteoff+join_wealth+join_gov+apol+cons+islamist+pray4+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey2)
plot(effect("BA_mil_3*Region", two, xlevels=list(coast_tun=c(0,1))),
     main="Willingness to fire on 2011 protesters (Tunisia)", 
     xlab="Ben Ali favored Sahel in promotions (5=str. agree)", 
     ylab="Would have fired (%)",
     colors="black", band.colors="grey0", alternating=F)


three <- lm((fire)~counterbalanced+gamal+eliteoff+join_wealth+join_gov+apol+cons+islamist+pray4+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural+registered, data=survey3)

plot_model(three, terms=c("counterbalanced","gamal","eliteoff","join_wealth","join_gov","cons"), type="est", title="Willingness to fire on 2011 protesters (Egypt)", show.p=T, show.values=T, p.threshold=c(0.10,0.05, 0.01), axis.labels=rev(c("Counterbalanced","Feared Gamal","From the Elite","Joined for Wealth","Joined for Power","Conscript")), value.offset=.3, dot.size=1, value.size=3.8, color="black") + ylim(-.2, .3) + set_theme(geom.outline.color = "antiquewhite4", base = theme_bw()) + geom_hline(yintercept=0, color="black", size=.2) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank()) + theme(axis.line = element_line(color = 'black')) +  theme(plot.title = element_text(hjust = 0.5))



##########
## Coup ##
##########

one <- lm(as.numeric(coup_ammar)~DO.BL.Prime_Tun+rebalanced+interior+eliteoff+join_wealth+join_gov+apol+cons+islamist+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey2)

two <- lm(as.numeric(coup_ammar)~DO.BL.Prime_Tun+rebalanced+affirm*interior+eliteoff+join_wealth+join_gov+apol+cons+islamist+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey2)

three <- lm(as.numeric(direct)~DO.BL.Prime+encroached+eliteoff+join_wealth+join_gov+apol+cons+islamist+shafik+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey3)

four <- lm(as.numeric(direct)~DO.BL.Prime+encroached+brotherhoodized+eliteoff+join_wealth+join_gov+apol+cons+islamist+shafik+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey3)

stargazer(one, two, three, four, single.row=T, digits=2, covariate.labels=c("Prime","Prime","Prime","Troika rebalanced","Troika elevated interior","from interior","Prime","Prime","Morsi encroached","Morsi Brotherhoodized","Elite Officer","Joined for Wealth","Joined for Power","Apolitical","Conscript","Islamist","Voted Shafik","Registered to Vote","Western Training","Joined to Defend Country","Joined for Thrill of Fighting","Deployed 2011","Active-Duty","Army","Navy","Soldier","NCO","Jun. Officer","Sen. Officer","Education","Rural","Troika elev. interior*from interior"), dep.var.labels=c("Troika (Tunisia)","Morsi (Egypt)"), dep.var.caption="Dependent Variable: Support for Coup against...", title="Support for a Military Coup in 2013")



## Plots 

two <- lm(as.numeric(coup_ammar)~DO.BL.Prime_Tun+rebalance+affirm*interior+eliteoff+join_wealth+join_gov+apol+cons+islamist+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey2)
plot(effect("rebalance", two), main="Support for a Coup in 2013 (Tunisia)", xlab="The Troika rebalanced the Military (5=str. agree)", ylab="Support for Coup (1-5)", colors="black", band.colors="grey0")

three <- lm(as.numeric(direct)~DO.BL.Prime+encroach+eliteoff+join_wealth+join_gov+apol+cons+islamist+shafik+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey3)
plot(effect("encroach", three), main="Support for the Coup in 2013 (Egypt)", xlab="Morsi encroached on military's interests (5=str. agree)", ylab="Support for Coup (0-1)", colors="black", band.colors="grey0")

four <- lm(as.numeric(direct)~DO.BL.Prime+encroach+morsi_mil_4+eliteoff+join_wealth+join_gov+apol+cons+islamist+shafik+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey3)
plot(effect("morsi_mil_4", four), main="Support for the Coup in 2013 (Egypt)", xlab="Morsi Brotherhoodized the military (5=str. agree)", ylab="Support for Coup (0-1)", colors="black", band.colors="grey0")





#################
## Priming Exp ##
#################


## Survey 1: Tunisia
data.summary <- data.frame(
  treatment=c("Worsened (N=37)", "Improved (N=35)"),
  mean=tapply(survey1$democracy, survey1$prime, mean, na.rm=TRUE),
  n=tapply(survey1$democracy, survey1$prime, length),
  sd=tapply(survey1$democracy, survey1$prime, sd, na.rm=TRUE))

data.summary$sem <- (data.summary$sd/sqrt(data.summary$n))
data.summary$me <- qt(1-.1/2, df=data.summary$n)*data.summary$sem

ggplot(data.summary, 
       aes(x=treatment, y=mean)) +  
  geom_bar(position=position_dodge(width=0.9), stat="identity", fill=c("gray40", "gray80")) + 
  geom_errorbar(aes(ymin=mean-me, ymax=mean+me, width=0.1),
                position=position_dodge(width=0.9)) +
  ggtitle("Priming Experiment") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(2.5, 4),breaks=seq(2.5, 4, by=0.5),oob=rescale_none) +
  xlab("") +
  ylab("Support Democracy") +
  geom_text(data=data.summary, 
            aes(x=treatment, y=mean,
            label=round(mean,2)),
            vjust=c(-4,-3),size=6,position=position_dodge(0.9)) +
  theme(text = element_text(size=17)) +
  geom_text(aes(x=1.5, y=3.9, label="p=0.056"), size=5)



## Survey 2: Tunisia
data.summary <- data.frame(
  treatment=c("Control (N=67)", "Ended Disc (N=68)","Improved (N=66)", "Worsened (N=70)"),
  mean=tapply(as.numeric(survey2$coup_ammar), survey2$DO.BL.Prime_Tun, mean, na.rm=TRUE),
  n=tapply(as.numeric(survey2$coup_ammar), survey2$DO.BL.Prime_Tun, length),
  sd=tapply(as.numeric(survey2$coup_ammar), survey2$DO.BL.Prime_Tun, sd, na.rm=TRUE))

data.summary$treatment <- factor(data.summary$treatment, levels=c("Control (N=67)", "Worsened (N=70)", "Improved (N=66)", "Ended Disc (N=68)"))
data.summary$sem <- (data.summary$sd/sqrt(data.summary$n))
data.summary$me <- qt(1-.1/2, df=data.summary$n)*data.summary$sem

ggplot(data.summary, 
       aes(x=treatment, y=mean)) +  
  geom_bar(position=position_dodge(width=0.9), stat="identity", fill=c("gray50", "gray40","gray50","gray80")) + 
  geom_errorbar(aes(ymin=mean-me, ymax=mean+me, width=0.1),
                position=position_dodge(width=0.9)) +
  ggtitle("Priming Experiment (Tunisia)") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(1, 3.5),breaks=seq(1, 3.5, by=1),oob=rescale_none) +
  xlab("") +
  ylab("Support for Coup (1-5)") +
  geom_text(data=data.summary, 
            aes(x=treatment, y=mean,
            label=round(mean,2)),
            vjust=c(-3.6),size=5,position=position_dodge(0.9)) +
  theme(text = element_text(size=15)) +
  geom_text(aes(x=1.5, y=3.4, label="p=0.096"), size=5)


## Survey 3: Egypt
data.summary <- data.frame(
  treatment=c("Control (N=88)", "Improved (N=109)", "Worsened (N=86)"),
  mean=tapply(as.numeric(survey3$direct[survey3$officer==1]), survey3$DO.BL.Prime[survey3$officer==1], mean, na.rm=TRUE),
  n=tapply(as.numeric(survey3$direct[survey3$officer==1]), survey3$DO.BL.Prime[survey3$officer==1], length),
  sd=tapply(as.numeric(survey3$direct[survey3$officer==1]), survey3$DO.BL.Prime[survey3$officer==1], sd, na.rm=TRUE))

data.summary$treatment <- factor(data.summary$treatment, levels=c("Control (N=88)", "Worsened (N=86)","Improved (N=109)"))
data.summary$sem <- (data.summary$sd/sqrt(data.summary$n))
data.summary$me <- qt(1-.1/2, df=data.summary$n)*data.summary$sem

ggplot(data.summary, 
       aes(x=treatment, y=mean)) +  
  geom_bar(position=position_dodge(width=0.9), stat="identity", fill=c("gray40", "gray50","gray80")) + 
  geom_errorbar(aes(ymin=mean-me, ymax=mean+me, width=0.1),
                position=position_dodge(width=0.9)) +
  ggtitle("Priming Experiment (Egypt-Officers)") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(limits=c(0, 1.1),breaks=seq(0, 1, by=0.2),oob=rescale_none) +
  xlab("") +
  ylab("Support for Coup (0-1)") +
  geom_text(data=data.summary, 
            aes(x=treatment, y=mean,
            label=round(mean,2)),
            vjust=c(-2.5),size=5,position=position_dodge(0.9)) +
  theme(text = element_text(size=16)) +
  geom_text(aes(x=1.5, y=1.06, label="p=0.02"), size=5)







###################
## Balance Plots ##
###################


## Survey 1: Tunisia

balance <- as.data.frame(rbind(
  summary(lm(rank3~prime, data=survey1))$coef[2,1:2],
  summary(lm(army~prime, data=survey1))$coef[2,1:2],
  summary(lm(fromSahel~prime, data=survey1))$coef[2,1:2],
  summary(lm(train_US~prime, data=survey1))$coef[2,1:2],
  summary(lm(islamist~prime, data=survey1))$coef[2,1:2]
))

balance$name <- c("Rank (1-3)", "Army (0-1)", "From Sahel (0-1)", "US-Trained (0-1)","Islamist (0-1)")
balance$name <- factor(balance$name, levels=c("Rank (1-3)", "Army (0-1)", "From Sahel (0-1)", "US-Trained (0-1)","Islamist (0-1)"))
balance$name <- factor(balance$name, levels = rev(levels(factor(balance$name))))

colnames(balance) <- c("mean", "se", "name")


ggplot(balance,
       aes(x=mean, y=name)) +
  geom_point(stat="identity", position="identity") +
  geom_errorbarh(aes(xmin=mean-1.96*se, 
                     xmax=mean+1.96*se, height=0.5), size=0.2) +
  geom_vline(aes(xintercept=0)) +
  ggtitle("Balance Plot (Priming Experiment)") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Treatment - Control") +
  ylab("") +
  scale_x_continuous(limits=c(-1, 1),breaks=seq(-1, 1, by=0.5)) +
  theme(text = element_text(size=17))


## Survey 2: Tunisia

balance <- as.data.frame(rbind(
  summary(lm(eliteoff~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(join_wealth~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(join_gov~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(apol~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(cons~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(islamist~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(registered~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(train_west~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(join_dut~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(join_adv~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(deploy~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(active~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm((branch=="Army")~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm((branch=="Navy")~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm((rank3=="Soldier")~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm((rank3=="NCO")~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm((rank3=="Jun. Officer")~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm((rank3=="Sen. Officer")~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(edu2~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2],
  summary(lm(rural~DO.BL.Prime_Tun, data=survey2))$coef[4,1:2]
))


balance$name <- c("From the Elite","Joined for Wealth","Joined for Power","Apolitical","Conscript","Islamist","Registered","Western Training","Joined to Defend Country","Joined for Thrill","Deployed 2011","Active-Duty","Army","Navy","Soldier","NCO","Jun. Officer","Sen. Officer","Education","Rural")
balance$name <- factor(balance$name, levels=c("From the Elite","Joined for Wealth","Joined for Power","Apolitical","Conscript","Islamist","Registered","Western Training","Joined to Defend Country","Joined for Thrill","Deployed 2011","Active-Duty","Army","Navy","Soldier","NCO","Jun. Officer","Sen. Officer","Education","Rural"))
balance$name <- factor(balance$name, levels = rev(levels(factor(balance$name))))

colnames(balance) <- c("mean", "se", "name")


ggplot(balance,
       aes(x=mean, y=name)) +
  geom_point(stat="identity", position="identity") +
  geom_errorbarh(aes(xmin=mean-1.96*se, 
                     xmax=mean+1.96*se, height=0.5), size=0.2) +
  geom_vline(aes(xintercept=0)) +
  ggtitle("Balance Plot (Tunisia)") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Worsened - Control") +
  ylab("") +
  scale_x_continuous(limits=c(-0.5, 0.5),breaks=seq(-0.5, 0.5, by=0.2)) +
  theme(text = element_text(size=15))




## Survey 3: Egypt

balance <- as.data.frame(rbind(
  summary(lm(eliteoff~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(join_wealth~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(join_gov~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(apol~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(cons~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(islamist~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(shafik~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(registered~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(train_west~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(join_dut~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(join_adv~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(deploy~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(active~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm((branch=="Army")~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm((branch=="Navy")~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm((rank3=="Soldier")~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm((rank3=="NCO")~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm((rank3=="Jun. Officer")~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm((rank3=="Sen. Officer")~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(edu2~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2],
  summary(lm(rural~DO.BL.Prime, data=survey3[survey3$officer==1,]))$coef[3,1:2]
))


balance$name <- c("From the Elite","Joined for Wealth","Joined for Power","Apolitical","Conscript","Islamist","Shafik","Registered","Western Training","Joined to Defend Country","Joined for Thrill","Deployed 2011","Active-Duty","Army","Navy","Soldier","NCO","Jun. Officer","Sen. Officer","Education","Rural")
balance$name <- factor(balance$name, levels=c("From the Elite","Joined for Wealth","Joined for Power","Apolitical","Conscript","Islamist","Shafik","Registered","Western Training","Joined to Defend Country","Joined for Thrill","Deployed 2011","Active-Duty","Army","Navy","Soldier","NCO","Jun. Officer","Sen. Officer","Education","Rural"))
balance$name <- factor(balance$name, levels = rev(levels(factor(balance$name))))
colnames(balance) <- c("mean", "se", "name")


ggplot(balance,
       aes(x=mean, y=name)) +
  geom_point(stat="identity", position="identity") +
  geom_errorbarh(aes(xmin=mean-1.96*se, 
                     xmax=mean+1.96*se, height=0.5), size=0.2) +
  geom_vline(aes(xintercept=0)) +
  ggtitle("Balance Plot (Egypt)") +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Worsened - Control") +
  ylab("") +
  scale_x_continuous(limits=c(-0.4, 0.4),breaks=seq(-0.4, 0.4, by=0.1)) +
  theme(text = element_text(size=15))





##############
## List exp ##
##############


one <- ictreg(list~encroached+eliteoff+join_wealth+join_gov+apol+cons+islamist+shafik+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural, data=survey3, treat="treatlist", J=4, method="lm")

two <- ictreg(list~encroached+eliteoff+join_wealth+join_gov+apol+cons+islamist+shafik+registered+train_west+join_dut+join_adv+deploy+active+branch+(rank3)+edu2+rural+brotherhoodized, data=survey3, treat="treatlist", J=4, method="lm")

one <- as.data.frame(cbind(one$par.treat, one$se.treat, 2*pnorm(-abs(one$par.treat/one$se.treat))))
one$p <- ifelse(one$V3>0.049 & one$V3<0.10, "*",
                ifelse(one$V3>0.01 & one$V3<0.05, "**",
                       ifelse(one$V3<0.01, "***", "")))
one$model1 <- sprintf("%2.2f%s (%2.2f)", one$V1, one$p, one$V2)
one <- one[5]

two <- as.data.frame(cbind(two$par.treat, two$se.treat, 2*pnorm(-abs(two$par.treat/two$se.treat))))
two$p <- ifelse(two$V3>0.049 & two$V3<0.10, "*",
                ifelse(two$V3>0.01 & two$V3<0.05, "**",
                       ifelse(two$V3<0.01, "***", "")))
two$model2 <- sprintf("%2.2f%s (%2.2f)", two$V1, two$p, two$V2)
two <- two[5]

data.summary <- cbind(rbind(one,c("")),two)
row.names(data.summary) <- c("Constant","Morsi encroached","Elite Officer","Joined for Wealth","Joined for Power","Apolitical","Conscript","Islamist","Shafik","Registered to Vote","Western Training","Joined to Defend Country","Joined for Thrill","Deployed 2011","Active-Duty","Army","Navy","Soldier","NCO","Jun. Officer","Sen. Officer","Education","Rural","Morsi Brotherhoodized")

xtable(data.summary, caption="Support for 2013 Coup in Egypt (List Experiment)")



